clear all; clc;close all;
%Updated: 05/07/2018

a   =load('agrid2.txt');
b   =load('bgrid2.txt');

vcon0=load('gvcon.txt');
vaut0=load('gvaut.txt');

p0  =load('gp.txt');
h0  =load('gh.txt');
ibp0=load('gibp.txt');
ct0 =load('gct.txt');
ca0 =load('gca.txt');
gn0 =load('ggn.txt');
tau0=load('gtau.txt');
w0  =load('gw.txt');
def0=load('gdef1.txt');
q0  =load('gq.txt');

paut0  =load('gpaut.txt');
haut0  =load('ghaut.txt');
ctaut0 =load('gctaut.txt');
gnaut0 =load('ggnaut.txt');
tauaut0=load('gtauaut.txt');
waut0  =load('gwaut.txt');

fm0     =load('fm_ss.txt');
rchange0=load('rchange_ss.txt');
derivb0=load('derivativeb_ss.txt');
bpvalue30=load('bpoptimal_ss.txt');
bpvalue20=load('bpgn2_ss.txt');
qvalue30=load('qoptimal_ss.txt');
qvalue20=load('qgn2_ss.txt');

% fm0      =load('fm_mean.txt');
% rchange0 =load('rchange_mean.txt');
% derivb0  =load('derivativeb_mean.txt');
% bpvalue30=load('bpoptimal_mean.txt');
% bpvalue20=load('bpgn2_mean.txt');
% qvalue30 =load('qoptimal_mean.txt');
% qvalue20 =load('qgn2_mean.txt');

na=rows(a);
nk=1;
ns=na/nk;
nb=rows(b);


alpha = 0.75;
pm     = 1;
gamma  = 0.43;
omegac = 0.3;
mu     = 1;
lambda = 0.184;
r      = 0.02;

atrans = ones(nb,1)*a';

for i=1:na;
    for j=1:nb;
        vcon(i,j)  =vcon0((i-1)*nb + j);
        
        p(i,j)  =p0((i-1)*nb + j);
        h(i,j)  =h0((i-1)*nb + j);
        ibp(i,j)=ibp0((i-1)*nb + j);
        ct(i,j) =ct0((i-1)*nb + j);
        ca(i,j) =ca0((i-1)*nb + j);
        gn(i,j) =gn0((i-1)*nb + j);       
        tau(i,j)=tau0((i-1)*nb + j);       
        w(i,j)  =w0((i-1)*nb + j);   
        def(i,j)=def0((i-1)*nb + j);   
        q(i,j)  =q0((i-1)*nb + j);   
        rec(i,j)=w(i,j)*h(i,j)*tau(i,j);
        gnr(i,j) = 0.200546+0.104336*a(i)-0.017938*b(j)/0.204;
        ratior(i,j) = 0.194311-0.06996*a(i)-0.03164*b(j)/0.204;
        ratio(i,j)=gn(i,j)/(h(i,j)^alpha);
        ratioby(i,j)=b(j)/exp(a(i));
        if (a(i)>-0.27982-1.21*b(j));
            gnfr(i,j) = 0.18357-0.29589*a(i)-0.256729*b(j);
%            gnfr(i,j) = 0.18357-0.31965*a(i)-0.261982*b(j)+0.294478*a(i)^2;
%            gnfr(i,j) = 0.181514-0.321359*a(i)-0.263489*b(j)+1.208048*a(i)^2-8.177788*a(i)^3;
            gnfr(i,j) = 0.1339-0.296282*a(i)-0.786244*b(j)+1.30132*a(i)^2-10.53735*a(i)^3-1.379999*b(j)^3;
        else
            gnfr(i,j) = 0.49144+0.786687*a(i)+1.09981*b(j);
            gnfr(i,j) = 0.525076+0.75233*a(i)+1.32279*b(j)+7.152719*a(i)^3-1.407359*b(j)^3;
        end;
        if (def(i,j)==0);ratio(i,j) =NaN;ratior(i,j) =NaN;gnr(i,j)=NaN;gnfr(i,j)=NaN;gn(i,j) =NaN;end;
            
    end;
end;


yn = h.^alpha;

for i=1:na-1;
    for j=1:nb;
        dgn(i,j) =(gn(i+1,j)-gn(i,j)<+0.001);
        dgnyn(i,j) =(gn(i+1,j)/yn(i+1,j)-gn(i,j)/yn(i,j)<+0.001);
        
        dh(i,j)  = (h(i,j)==1);
    end;
end;

for i=1:na-1
    j=1;
    value=0;
while (value==0);
    j=j+1;
    ithres(i)=j;
    value = dgn(i,j);
    %disp([i,j,value,ithres(i)]);
end;
end;

thres_vec = [a(2:na),ithres'];

bp = b(ibp);


for i=1:na;
    vaut(i)  =vaut0(i);
    
    paut(i)  =paut0(i);
    haut(i)  =haut0(i);
    ctaut(i) =ctaut0(i);
    gnaut(i) =gnaut0(i);       
    tauaut(i)=tauaut0(i);       
    waut(i)  =waut0(i);   
    recaut(i)=waut(i)*haut(i)*tauaut(i);
	wwaut(i) = alpha*paut(i).*haut(i).^(alpha-1); 
    
end;

cn       = h.^alpha-gn;
cnaut    = haut.^alpha-gnaut;
const    = (omegac.*ct.^(-mu)+(1-omegac).*cn.^(-mu)).^(-1/mu);
constaut = (omegac.*ctaut.^(-mu)+(1-omegac).*cnaut.^(-mu)).^(-1/mu);

for i=1:na;
    temp1 = find(ct(i,:) > 0);
    temp2 = find(cn(i,:) > 0);
    temp3 = find(p(i,:) > 0);

    iminct(i) = temp1(1);
    imincn(i) = temp2(1);
    iminp(i) = temp3(1);
    
    imin(i)  = max(iminct(i),max(imincn(i),iminp(i)));
    
    v(i,1:imin(i)-1)  = NaN; 
    p(i,1:imin(i)-1)  = NaN; 
    bp(i,1:imin(i)-1) = NaN; 
    ibp(i,1:imin(i)-1)= NaN; 
    ctE(i,1:imin(i)-1) = NaN; 
    ctU(i,1:imin(i)-1) = NaN; 
    cnE(i,1:imin(i)-1) = NaN; 
    cnU(i,1:imin(i)-1) = NaN; 
    constE(i,1:imin(i)-1) = NaN;     
    constU(i,1:imin(i)-1) = NaN;     
    gn(i,1:imin(i)-1) = NaN; 
    h(i,1:imin(i)-1)  = NaN; 
    ca(i,1:imin(i)-1)  = NaN; 
    tau(i,1:imin(i)-1)= NaN; 
	rec(i,1:imin(i)-1)= NaN;
	ww(i,1:imin(i)-1)= NaN;
	q(i,1:imin(i)-1)  = NaN;
  
end;


for i=1:ns;
for j=1:nk;        
    vcon2(i,j,:)  = vcon((i-1)*nk+j,:); 
    p2(i,j,:)  = p((i-1)*nk+j,:);  
    bp2(i,j,:) = bp((i-1)*nk+j,:);  
    ibp2(i,j,:)= ibp((i-1)*nk+j,:);  
    ct2(i,j,:) = ct((i-1)*nk+j,:);  
    cn2(i,j,:) = cn((i-1)*nk+j,:);  
    const2(i,j,:) = const((i-1)*nk+j,:);         
    gn2(i,j,:) = gn((i-1)*nk+j,:);  
	ca2(i,j,:) = ca((i-1)*nk+j,:);  
    h2(i,j,:)  = h((i-1)*nk+j,:);   
    tau2(i,j,:)= tau((i-1)*nk+j,:);  
	rec2(i,j,:)= rec((i-1)*nk+j,:); 
    def2(i,j,:)= def((i-1)*nk+j,:);  
	q2(i,j,:)  = q((i-1)*nk+j,:); 
	ww2(i,j,:) = alpha*p2(i,j,:).*h2(i,j,:).^(alpha-1); 
    

    vaut2(i,j)  = vaut((i-1)*nk+j); 
    paut2(i,j)  = paut((i-1)*nk+j);  
    ctaut2(i,j) = ctaut((i-1)*nk+j);  
    cnaut2(i,j) = cnaut((i-1)*nk+j);  
    constaut2(i,j) = constaut((i-1)*nk+j);          
    gnaut2(i,j) = gnaut((i-1)*nk+j);  
    haut2(i,j)  = haut((i-1)*nk+j);  
    tauaut2(i,j)= tauaut((i-1)*nk+j);  
	wwaut2(i,j) = alpha*paut2(i,j).*haut2(i,j).^(alpha-1); 
end;
end;

nettax2    = p2.*gn2-ca2;           %net taxes
nettaxaut2 = paut2.*gnaut2;         %net taxes

fiscal2 = tau2 - p2.*gn2;           %fiscal surplus
fiscalaut2=tauaut2-paut2.*gnaut2;   %fiscal surplus

nfiscal2 = nettax2 - p2.*gn2;           %fiscal surplus net of transfers
nfiscalaut2=nettaxaut2-paut2.*gnaut2;   %fiscal surplus net of transfers

lowa = round(2*ns/6);
meda = round(ns/2);
hgha = round(4*ns/6);

b =-b/(lambda+r);  %change of sign
a = exp(a);

% %Compute fiscal multipliers an borrowing costs
% for i=1:nb;
% for j=1:na;
% 	fm(j,i)      =fm0((i-1)*na + j);     
%   	rchange(j,i) =rchange0((i-1)*na + j);     
%   	derivb(j,i)  =-derivb0((i-1)*na + j);     
%   	bpvalue2(j,i)=bpvalue20((i-1)*na + j);     
%   	bpvalue3(j,i)=bpvalue30((i-1)*na + j);     
%   	qvalue2(j,i)=qvalue20((i-1)*na + j);     
%   	qvalue3(j,i)=qvalue30((i-1)*na + j);     
%     if (def(j,i)==0);fm(j,i)=NaN;rchange(j,i)=NaN;derivb(j,i)=NaN;bpvalue2(j,i)=NaN;...
%             bpvalue3(j,i)=NaN;qvalue2(j,i)=NaN;qvalue3(j,i)=NaN;end;     
% end;
% end;
% 
% %debt indices
% %default IR: 
% %%=122 mean debt, =196 0.5*mean debt,   =159 0.75*mean debt, 
% 
% %risk-free debt IR:
% %%=1888 mean debt, =1944 0.5*mean debt, =1916 0.75*mean debt,  =1774 2*meandebt 
% %=1973 0.25*mean debt, =2001 zero,
% 
% lowb = 159;     %=171 0.5*mean debt, =116 0.75*mean debt
% hghb = 122;      %mean debt
% 
% figure;
% subplot(2,1,1);
% plot(b,fm(lowa,:),'--r',b,fm(hgha,:),'-b','LineWidth',1.5);
% legend('low y^T','high y^T','Location','NorthWest','Orientation','vertical');xlabel('debt');
% AX1 = gca;box on;
% set(AX1,'YTick',-1:1:3);
% axis(AX1, [-0.5 2 -1 3]);
% title('fiscal mutipliers');xlabel('debt');
% subplot(2,1,2);
% plot(b,derivb(lowa,:),'--r',b,derivb(hgha,:),'-b','LineWidth',1.5);
% legend('low y^T','high y^T','Location','NorthWest','Orientation','vertical');xlabel('debt');%ylabel('%');
% title('borrowing costs');xlabel('debt');
% AX1 = gca;box on;
% set(AX1,'YTick',0:0.5:1.5);
% axis(AX1, [-0.5 2 -0.1 1.5]);
% 
% figure;
% subplot(2,1,1);
% plot(a,fm(:,lowb),'--r',a,fm(:,hghb),'-b','LineWidth',1.5);
% legend('low debt','high debt','Location','NorthEast','Orientation','vertical');xlabel('y^T');
% AX1 = gca;box on;
% set(AX1,'YTick',-1:1:3);
% axis(AX1, [0.85 1.15 -1 3]);
% title('fiscal mutipliers');xlabel('debt');
% subplot(2,1,2);
% plot(a,derivb(:,lowb),'--r',a,derivb(:,hghb),'-b','LineWidth',1.5);
% legend('low debt','high debt','Location','NorthEast','Orientation','vertical');xlabel('y^T');%ylabel('%');
% title('borrowing costs');xlabel('debt');
% AX1 = gca;box on;
% set(AX1,'YTick',0:0.5:1);
% axis(AX1, [0.85 1.15 -0.1 1]);

[i1,i2]=find(def(lowa,:)>0.5);
start=i2(1);

figure;
subplot(2,1,1);
plot(b(start:end),gn(lowa,start:end)-gn(meda,start:end),'-b','LineWidth',1.5);xlabel('debt');title('Response on impact of g^N of y^T adverse shock');
subplot(2,1,2);
plot(b(start:end),bp(lowa,start:end)-bp(meda,start:end),'-b','LineWidth',1.5);xlabel('debt');title('Response on impact of debt of y^T adverse shock');

%put together policy fcns of repayment and autarky
mat_ones = ones(size(def));
risk_free_debt=1;

if (risk_free_debt==1);
    p3  =p;
    h3  =h;
    ct3=ct;
    cn3=cn;
    ca3 =ca;
    gn3 =gn;       
    tau3=tau;       
    w3  =w;   
    rec3=rec;
else;
    p3  =def.*p+(mat_ones-def).*(paut'*ones(1,nb));
    h3  =def.*h+(mat_ones-def).*(haut'*ones(1,nb));
    ct3=def.*ct+(mat_ones-def).*(ctaut'*ones(1,nb));
    cn3=def.*cn+(mat_ones-def).*(cnaut'*ones(1,nb));
    ca3 =def.*ca+(mat_ones-def).*0;
    gn3 =def.*gn+(mat_ones-def).*(gnaut'*ones(1,nb));       
    tau3=def.*tau+(mat_ones-def).*(tauaut'*ones(1,nb));       
    w3  =def.*w+(mat_ones-def).*(waut'*ones(1,nb));   
    rec3=def.*rec+(mat_ones-def).*(recaut'*ones(1,nb));
end;

%cnaut3 = haut3.*cnEaut3 + (1-haut3).*cnUaut3;
%ctaut3 = haut3.*ctEaut3 + (1-haut3).*ctUaut3;

[i1,i2] = find(def>0);
q3  =NaN; bp3 =NaN; sp3 =NaN; gn3 = NaN; h3=NaN; p3=NaN;

for i=1:na;
for j=1:nb;
    if (def(i,j)>0); q3(i,j)=q(i,j);bp3(i,j)=bp(i,j);sp3(i,j)=100*((1+q3(i,j).^(-1)-lambda)./(1+r)-1);p3=p;h3=h;gn3=gn;tau3=tau;w3=w;rec3=rec;
    else;
        q3(i,j)=NaN;bp3(i,j)=NaN;sp3(i,j)=NaN;p3=NaN;h3=NaN;gn3=NaN;tau3=NaN;w3=NaN;rec3=NaN;
    end;
end;
end;

nettax3    = p3.*gn3-ca3;           %net taxes
fiscal3 = tau3 - p3.*gn3;           %fiscal surplus
nfiscal3 = nettax3 - p3.*gn3;           %fiscal surplus net of transfers

%default decision and prices
def21(:,:) = def2(:,1,:);
def22(:,:) = def2(:,1,:);

price21(:,:) = q2(:,1,:);
price22(:,:) = q2(:,1,:);


% figure;
% plot(b(1:nb),ctE21(lowa,:),'-r',b(1:nb),ctE21(meda,:),':r',b(1:nb),ctE21(hgha,:),'--r',...
%     b(1:nb),ctU21(lowa,:),'-b',b(1:nb),ctU21(meda,:),':b',b(1:nb),ctU21(hgha,:),'--b','LineWidth',1.5);
% xlabel('b');legend('c^T_E y^T low','c^T_E y^T med','c^T_E y^T hgh','c^T_U y^T low','c^T_U y^T med','c^T_U y^T hgh');
% xlim([min(b) 1]);
% 
% figure;
% subplot(2,2,1);
% contourf(b(1:nb),a(1:ns),def21(:,:),20);
% xlabel('b');
% ylabel('a');
% title('default decision \kappa_{L}');
% subplot(2,2,2);
% contourf(b(1:nb),a(1:ns),def22(:,:),20);
% xlabel('b');
% ylabel('a');
% title('default decision \kappa_{H}');
% subplot(2,2,[3 4]);
% plot(b(1:nb),price21(lowa,:),'-r',b(1:nb),price21(meda,:),':r',b(1:nb),price21(hgha,:),'--r',...
%     b(1:nb),price22(lowa,:),'-b',b(1:nb),price22(meda,:),':b',b(1:nb),price22(hgha,:),'--b','LineWidth',1.5);
% xlabel('b');
% 
% 
% %%
% %value functions: repayment vs autarky
% vcon21(:,:) = vcon2(:,1,:);
% vcon22(:,:) = vcon2(:,1,:);
% vaut21(:,:) = vaut2(:,1);
% vaut22(:,:) = vaut2(:,1);
% 
% figure;
% plot(b(1:nb),vcon21(lowa,:),'-r',b(1:nb),vcon21(meda,:),':r',b(1:nb),...
%     vaut21(lowa,:)*ones(size(b)),'-b',b(1:nb),vaut21(meda)*ones(size(b)),':b','LineWidth',1.5);
% xlabel('b');title('value functions');
% ylim([-16 -13]);legend('rep low y^T','rep high y^T','aut low y^T','aut high y^T')

%% Plot variables as function of yT for two debt levels
%optimal policies and prices: repayment vs autarky
%NameArray = {'LineStyle','Color','LineWidth'};
%ValueArray = {'-',':','-',':';'red','red','blue','blue';1.5,1.5,1.5,1.5}';

%NameArray0 = {'LineStyle','Color','LineWidth'};
%ValueArray0 = {'-','-';'red','blue';1.5,1.5}';
NameArray = {'LineStyle','Color','LineWidth'};
NameArray0 = {'LineStyle','Color','LineWidth'};
ValueArray0 = {'--','-';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';
ValueArray = {'--','-';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';
ValueArray2 = {'--','-','--','-';[1  0.4  0.4],[0.078  0.1686  0.549],[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';
ValueArray4 = {'--','--','-','-';'red','blue','red','blue';1.5,1.5,1.5,1.5}';
ValueArray5 = {':',':';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';

endb = 2;
lowb = 159;
hghb = 122;

figure;
subplot(5,2,1);P = plot(a,100*squeeze(1-h3(:,lowb)),a,100*squeeze(1-h3(:,hghb)));
title('unemployment','Fontweight','normal');
AX1 = gca;box on;ylabel('%');
set(AX1,'YTick',0:20:40);
axis(AX1, [a(1) a(end) -2 40]);
set(P,NameArray,ValueArray);
subplot(5,2,2);P = plot(a,squeeze(p3(:,lowb)),a,squeeze(p3(:,hghb)));
title('p^{N}','Fontweight','normal'); 
%lgd=legend('low y^T','high y^T');lgd.Location='SouthEast';
legend('low debt','high debt','Location','SouthEast','Orientation','vertical');
AX1 = gca;box on;
set(AX1,'YTick',3.5:0.5:4.5);
axis(AX1, [a(1) a(end) 3.5 4.5]);
set(P,NameArray,ValueArray);
subplot(5,2,3);P = plot(a,squeeze(gn3(:,lowb)),a,squeeze(gn3(:,hghb)));
title('g^{N}','Fontweight','normal'); 
AX1 = gca;box on;
set(AX1,'YTick',0.1:0.1:0.3);
axis(AX1, [a(1) a(end) 0.1 0.3]);
set(P,NameArray,ValueArray);
subplot(5,2,4);P = plot(a,-squeeze(bp3(:,lowb))/(lambda+r),a,-squeeze(bp3(:,hghb))/(lambda+r));
title('debt','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.8:0.2:1.2);
axis(AX1, [a(1) a(end) 0.8 1.2]);
set(P,NameArray0,ValueArray0);
subplot(5,2,5);P = plot(a,squeeze(cn3(:,lowb)),a,squeeze(cn3(:,hghb)));
title('c^{N}','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.5:0.2:0.9);
axis(AX1, [a(1) a(end) 0.5 0.9]);
set(P,NameArray,ValueArray);
subplot(5,2,6);P = plot(a,squeeze(ct3(:,lowb)),a,squeeze(ct3(:,hghb)));
title('c^{T}','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.7:0.2:1.1);
axis(AX1, [a(1) a(end) 0.7 1.1]);
set(P,NameArray,ValueArray);
subplot(5,2,7);P = plot(a,squeeze(w3(:,lowb)),a,squeeze(w3(:,hghb)));
title('wages','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',2.5:0.5:3.5);
axis(AX1, [a(1) a(end) 2.5 3.5]);
set(P,NameArray,ValueArray);
%version 1: with bond spreads
subplot(5,2,8);P = plot(a,squeeze(sp3(:,lowb)),a,squeeze(sp3(:,hghb)));
title('spreads','Fontweight','normal');
set(P,NameArray0,ValueArray0);ylabel('%');
AX1 = gca;box on;
set(AX1,'YTick',0:4:8);
axis(AX1, [a(1) a(end) -0.5 8]);
%version 2: with current account
subplot(5,2,9);P = plot(a,-squeeze(ca3(:,lowb)),a,-squeeze(ca3(:,hghb)));
%xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
title('CA balance','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',-0.1:0.1:0.2);
axis(AX1, [a(1) a(end) -0.1 0.2]);
set(P,NameArray0,ValueArray0);
%version 3: with (net) taxes
subplot(5,2,10);P = plot(a,squeeze(tau3(:,lowb)),a,squeeze(tau3(:,hghb)));
set(P,NameArray,ValueArray);
hold on;
PP = plot(a,squeeze(nettax3(:,lowb)),a,squeeze(nettax3(:,hghb)));
title('(net) taxes','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.6:0.3:1.2);
axis(AX1, [a(1) a(end) 0.6 1.2]);
set(PP,NameArray,ValueArray5);
%set(P,NameArray0,ValueArray0);
%subplot(6,2,11);P = plot(b,squeeze(fiscal2(lowa,1,:)),b,squeeze(fiscala2(lowa,1,:)),b,squeeze(fiscalr2(hgha,1,:)),b,squeeze(fiscala2(hgha,1,:)));
%title('primary surplus'); %xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
%subplot(4,2,8);P = plot(b,squeeze(taur2(lowa,1,:)),b,squeeze(taua2(lowa,1,:)),b,squeeze(taur2(hgha,1,:)),b,squeeze(taua2(hgha,1,:)));
%title('T'); xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
%set(P,NameArray,ValueArray);

dbstop;

%% Plot variables as function of debt for two yT levels
hr2  = h2;
cnr2 = cn2;
gnr2 = gn2;
ctr2 = ct2;
pr2  = p2;
taur2= tau2;
nettaxr2= nettax2;
fiscalr2  = fiscal2;
nfiscalr2 = nfiscal2;
bpr2 = bp2;
wwr2 = ww2;
qr2  = q2;
spr2 = 100*((1+qr2.^(-1)-lambda)./(1+r)-1);
car2 = ca2;

%optimal policies and prices: repayment vs autarky
%NameArray = {'LineStyle','Color','LineWidth'};
%ValueArray = {'-',':','-',':';'red','red','blue','blue';1.5,1.5,1.5,1.5}';

%NameArray0 = {'LineStyle','Color','LineWidth'};
%ValueArray0 = {'-','-';'red','blue';1.5,1.5}';
NameArray = {'LineStyle','Color','LineWidth'};
ValueArray = {'--','--','-','-';[1  0.4  0.4],[1  0.4  0.4],[0.078  0.1686  0.549],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';
ValueArray2 = {'--','-','--','-';[1  0.4  0.4],[0.078  0.1686  0.549],[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';
ValueArray4 = {'--','--','-','-';'red','blue','red','blue';1.5,1.5,1.5,1.5}';
%ValueArray5 = {':',':',':',':';'red','blue','red','blue';1.5,1.5,1.5,1.5}';
%ValueArray5 = {':',':',':',':';[1  0.4  0.4],[0.078  0.1686  0.549],[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';
ValueArray5 = {':',':',':',':';[1  0.4  0.4],[1  0.4  0.4],[0.078  0.1686  0.549],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';

NameArray0 = {'LineStyle','Color','LineWidth'};
ValueArray0 = {'--','-';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';

[ilowa_aut]=find(def2(lowa,1,:)==0);
[ilowa_rep]=find(def2(lowa,1,:)==1);
hr2(lowa,1,ilowa_aut) = NaN;
ha2(lowa,1,1:nb) = haut2(lowa,1);
ha2(lowa,1,ilowa_rep) = NaN;

[ihgha_aut]=find(def2(hgha,1,:)==0);
[ihgha_rep]=find(def2(hgha,1,:)==1);
hr2(hgha,1,ihgha_aut) = NaN;
ha2(hgha,1,1:nb) = haut2(hgha,1);
ha2(hgha,1,ihgha_rep) = NaN;

pr2(lowa,1,ilowa_aut) = NaN;
pa2(lowa,1,1:nb) = paut2(lowa,1);
pa2(lowa,1,ilowa_rep) = NaN;

pr2(hgha,1,ihgha_aut) = NaN;
pa2(hgha,1,1:nb) = paut2(hgha,1);
pa2(hgha,1,ihgha_rep) = NaN;

gnr2(lowa,1,ilowa_aut) = NaN;
gna2(lowa,1,1:nb) = gnaut2(lowa,1);
gna2(lowa,1,ilowa_rep) = NaN;

gnr2(hgha,1,ihgha_aut) = NaN;
gna2(hgha,1,1:nb) = gnaut2(hgha,1);
gna2(hgha,1,ihgha_rep) = NaN;

bpr2(lowa,1,ilowa_aut) = NaN;
bpr2(hgha,1,ihgha_aut) = NaN;

cnr2(lowa,1,ilowa_aut) = NaN;
cna2(lowa,1,1:nb) = cnaut2(lowa,1);
cna2(lowa,1,ilowa_rep) = NaN;

cnr2(hgha,1,ihgha_aut) = NaN;
cna2(hgha,1,1:nb) = cnaut2(hgha,1);
cna2(hgha,1,ihgha_rep) = NaN;

ctr2(lowa,1,ilowa_aut) = NaN;
cta2(lowa,1,1:nb) = ctaut2(lowa,1);
cta2(lowa,1,ilowa_rep) = NaN;

ctr2(hgha,1,ihgha_aut) = NaN;
cta2(hgha,1,1:nb) = ctaut2(hgha,1);
cta2(hgha,1,ihgha_rep) = NaN;

wwr2(lowa,1,ilowa_aut) = NaN;
wwa2(lowa,1,1:nb) = wwaut2(lowa,1);
wwa2(lowa,1,ilowa_rep) = NaN;

wwr2(hgha,1,ihgha_aut) = NaN;
wwa2(hgha,1,1:nb) = wwaut2(hgha,1);
wwa2(hgha,1,ihgha_rep) = NaN;

taur2(lowa,1,ilowa_aut) = NaN;
taua2(lowa,1,1:nb) = tauaut2(lowa,1);
taua2(lowa,1,ilowa_rep) = NaN;

taur2(hgha,1,ihgha_aut) = NaN;
taua2(hgha,1,1:nb) = tauaut2(hgha,1);
taua2(hgha,1,ihgha_rep) = NaN;

nettaxr2(lowa,1,ilowa_aut) = NaN;
nettaxa2(lowa,1,1:nb) = nettaxaut2(lowa,1);
nettaxa2(lowa,1,ilowa_rep) = NaN;

nettaxr2(hgha,1,ihgha_aut) = NaN;
nettaxa2(hgha,1,1:nb) = nettaxaut2(hgha,1);
nettaxa2(hgha,1,ihgha_rep) = NaN;

fiscalr2(lowa,1,ilowa_aut) = NaN;
fiscala2(lowa,1,1:nb) = fiscalaut2(lowa,1);
fiscala2(lowa,1,ilowa_rep) = NaN;

fiscalr2(hgha,1,ihgha_aut) = NaN;
fiscala2(hgha,1,1:nb) = fiscalaut2(hgha,1);
fiscala2(hgha,1,ihgha_rep) = NaN;

nfiscalr2(lowa,1,ilowa_aut) = NaN;
nfiscala2(lowa,1,1:nb) = nfiscalaut2(lowa,1);
nfiscala2(lowa,1,ilowa_rep) = NaN;

nfiscalr2(hgha,1,ihgha_aut) = NaN;
nfiscala2(hgha,1,1:nb) = nfiscalaut2(hgha,1);
nfiscala2(hgha,1,ihgha_rep) = NaN;

car2(lowa,1,ilowa_aut) = NaN;
car2(hgha,1,ihgha_aut) = NaN;

qr2(lowa,1,ilowa_aut) = NaN;
qr2(hgha,1,ihgha_aut) = NaN;

spr2(lowa,1,ilowa_aut) = NaN;
spr2(hgha,1,ihgha_aut) = NaN;

endb = 2;

figure;
subplot(5,2,1);P = plot(b,100*squeeze(1-hr2(lowa,1,:)),b,100*squeeze(1-ha2(lowa,1,:)),b,100*squeeze(1-hr2(hgha,1,:)),b,100*squeeze(1-ha2(hgha,1,:)));
title('unemployment','Fontweight','normal');
AX1 = gca;box on;ylabel('%');
set(AX1,'YTick',0:10:30);
axis(AX1, [min(b) endb -1 30]);
set(P,NameArray,ValueArray);
subplot(5,2,2);P = plot(b,squeeze(pr2(lowa,1,:)),b,squeeze(pr2(hgha,1,:)),b,squeeze(pa2(lowa,1,:)),b,squeeze(pa2(hgha,1,:)));
title('p^{N}','Fontweight','normal'); 
%lgd=legend('low y^T','high y^T');lgd.Location='SouthEast';
legend('low y^T','high y^T','Location','NorthEast','Orientation','vertical');
AX1 = gca;box on;
set(AX1,'YTick',3.5:0.5:5.5);
axis(AX1, [min(b) endb 3.5 5.5]);
set(P,NameArray,ValueArray2);
subplot(5,2,3);P = plot(b,squeeze(gnr2(lowa,1,:)),b,squeeze(gna2(lowa,1,:)),b,squeeze(gnr2(hgha,1,:)),b,squeeze(gna2(hgha,1,:)));
title('g^{N}','Fontweight','normal'); 
AX1 = gca;box on;
set(AX1,'YTick',0.12:0.08:0.28);
axis(AX1, [min(b) endb 0.12 0.28]);
set(P,NameArray,ValueArray);
%subplot(5,2,4);P = plot(b,squeeze(bpr2(lowa,1,:)),b,squeeze(bpr2(hgha,1,:)));
subplot(5,2,4);P = plot(b,-squeeze(bpr2(lowa,1,:))/(lambda+r),b,-squeeze(bpr2(hgha,1,:))/(lambda+r));
title('debt','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',-1.6:0.8:1.6);
axis(AX1, [min(b) endb -1.6 1.6]);
set(P,NameArray0,ValueArray0);
subplot(5,2,5);P = plot(b,squeeze(cnr2(lowa,1,:)),b,squeeze(cna2(lowa,1,:)),b,squeeze(cnr2(hgha,1,:)),b,squeeze(cna2(hgha,1,:)));
title('c^{N}','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.6:0.1:0.9);
axis(AX1, [min(b) endb 0.6 0.9]);
set(P,NameArray,ValueArray);
subplot(5,2,6);P = plot(b,squeeze(ctr2(lowa,1,:)),b,squeeze(cta2(lowa,1,:)),b,squeeze(ctr2(hgha,1,:)),b,squeeze(cta2(hgha,1,:)));
title('c^{T}','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.8:0.2:1.4);
axis(AX1, [min(b) endb 0.8 1.4]);
set(P,NameArray,ValueArray);
subplot(5,2,7);P = plot(b,squeeze(wwr2(lowa,1,:)),b,squeeze(wwa2(lowa,1,:)),b,squeeze(wwr2(hgha,1,:)),b,squeeze(wwa2(hgha,1,:)));
title('wages','Fontweight','normal');xlim([min(b) max(b)]);
AX1 = gca;box on;
set(AX1,'YTick',3:0.5:4.5);
axis(AX1, [min(b) endb 3 4.5]);
set(P,NameArray,ValueArray);
%version 1: with bond spreads
subplot(5,2,8);P = plot(b,squeeze(spr2(lowa,1,:)),b,squeeze(spr2(hgha,1,:)));
title('spreads','Fontweight','normal');
set(P,NameArray0,ValueArray0);ylabel('%');
AX1 = gca;box on;
set(AX1,'YTick',0:20:40);
axis(AX1, [min(b) endb -2 40]);
%version 2: with current account
subplot(5,2,9);P = plot(b,-squeeze(car2(lowa,1,:)),b,-squeeze(car2(hgha,1,:)));
%xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
title('CA balance','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',-0.3:0.3:0.3);
axis(AX1, [min(b) endb -0.3 0.3]);
set(P,NameArray0,ValueArray0);
%version 3: with (net) taxes
subplot(5,2,10);P = plot(b,squeeze(taur2(lowa,1,:)),b,squeeze(taua2(lowa,1,:)),b,squeeze(taur2(hgha,1,:)),b,squeeze(taua2(hgha,1,:)));
set(P,NameArray,ValueArray);
hold on;
PP = plot(b,squeeze(nettaxr2(lowa,1,:)),b,squeeze(nettaxa2(lowa,1,:)),b,squeeze(nettaxr2(hgha,1,:)),b,squeeze(nettaxa2(hgha,1,:)));
%xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
title('(net) taxes','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.25:0.5:1.25);
axis(AX1, [min(b) endb 0.25 1.25]);
set(PP,NameArray,ValueArray5);
%set(P,NameArray0,ValueArray0);
%subplot(6,2,11);P = plot(b,squeeze(fiscal2(lowa,1,:)),b,squeeze(fiscala2(lowa,1,:)),b,squeeze(fiscalr2(hgha,1,:)),b,squeeze(fiscala2(hgha,1,:)));
%title('primary surplus'); %xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
%subplot(4,2,8);P = plot(b,squeeze(taur2(lowa,1,:)),b,squeeze(taua2(lowa,1,:)),b,squeeze(taur2(hgha,1,:)),b,squeeze(taua2(hgha,1,:)));
%title('T'); xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
%set(P,NameArray,ValueArray);

dbstop;


%% plot time series
seriesp   = load('seriesp.txt');
seriesr   = load('seriesr.txt');
seriesh   = load('seriesh.txt');
seriesb   = load('seriesb.txt');
seriesct  = load('seriesct.txt');
seriesm   = 1; %load('seriesm.txt');
seriesgn  = load('seriesgn.txt');
seriestau = load('seriestau.txt');    
seriesa   = load('seriesa.txt');    
seriesbor = load('seriesbor.txt');    

n = length(seriesp);

seriescn= seriesh.^alpha-seriesgn;
seriesy = seriesa.*seriesm.^gamma+seriesp.*seriesh.^alpha; % -pm*seriesm;
seriesgny = seriesgn./seriesy;

time =1:n;

ValueArray1 = {'-';'blue';1.5}';
ValueArray2 = {'-','--';'red','blue';1.5,1.5}';
ValueArray3 = {'-','--','-';'red','blue','black';1.5,1.5,1.5}';

% No shaded areas for default episodes
% figure;
% subplot(3,2,1); P=plot(time,seriesa');title('a');
% set(P,NameArray,ValueArray1);
% subplot(3,2,2);P=plot(time,seriesp);title('p');
% set(P,NameArray,ValueArray1);
% subplot(3,2,3);P=plot(time,seriesct);
% %AX=legend('c_{T}','Location','NorthWest','Orientation','Horizontal');
% title('c_T');
% set(P,NameArray,ValueArray1);
% %LEG = findobj(AX,'type','text');
% %set(LEG,'FontSize',7);
% subplot(3,2,4); PPP=plot(time,seriescn,time,seriesh,time,seriesgn);
% axis([1 n -0.05 1.05]);
% AX=legend('h','c_{N}','g_{N}','Location','NorthWest','Orientation','Horizontal');title('h,c_{N} & g_{N}');
% set(PPP,NameArray,ValueArray3);
% LEG = findobj(AX,'type','text');
% set(LEG,'FontSize',7);
% subplot(3,2,5); P=plot(time,seriesr);title('spreads');
% set(P,NameArray,ValueArray1);axis([1 n -0.01 0.1]);
% subplot(3,2,6); P=plot(time,-seriesb);title('b');
% set(P,NameArray,ValueArray1);


seriesbor = 1-seriesbor;
seriesbor(199)=1;
seriesbor(200)=0;
seriesh(200) = 0;
seriesr(200) = 0.1;
ntime = 200;

figure;
subplot(3,2,1); 
shadedTimeSeries(time,seriesa',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('a');
AX1 = gca;box on;
set(AX1,'YTick',0.85:0.1:1.15);
axis(AX1, [time(1) time(ntime) 0.85 1.15]);
subplot(3,2,2);
shadedTimeSeries(time,seriesp',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('p');
AX1 = gca;box on;
set(AX1,'YTick',0:0.2:1);
axis(AX1, [time(1) time(ntime) 0.95*min(seriesp) 1.05*max(seriesp)]);
subplot(3,2,3);
shadedTimeSeries(time,seriesct',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('c_T');
AX1 = gca;box on;
set(AX1,'YTick',0.05:0.05:1.15);
axis(AX1, [time(1) time(ntime) 0.9*min(seriesct)  1.1*max(seriesct)]); 
subplot(3,2,4); 
shadedTimeSeries(time,seriesh',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('h,c_N & g_N');
AX1 = gca;box on;
axis(AX1, [time(1) time(ntime) 0 1.05]);
hold on;
H0 = line(time,seriesh','Color','b','LineWidth',2);
H1 = line(time,seriescn','Color','r','LineWidth',2,'LineStyle','--');
H2 = line(time,seriesgn','Color','g','LineWidth',2,'LineStyle',':');
hh = [H0, H1,H2];
legend(hh,'h','c_N','g_N','Location','SouthEast','Orientation','Horizontal');
AX1 = gca;box on;
axis(AX1, [time(1) time(ntime) 0 1.05]);
subplot(3,2,5); 
shadedTimeSeries(time,seriesr',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('spreads');
AX1 = gca;box on;
set(AX1,'YTick',0.0:0.025:0.1);
axis(AX1, [time(1) time(ntime) 0 0.1]);
subplot(3,2,6); P=plot(time,-seriesb);title('b');
shadedTimeSeries(time,-seriesb',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('b');
AX1 = gca;box on;
%set(AX1,'YTick',0.0:0.02:ceil(max(-1.1*seriesb/0.001))*0.001);
axis(AX1, [time(1) time(ntime) 0 max(-seriesb)]);

dbstop;